1 Loading in Packages + Dataset

## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpQoXYeb/downloaded_packages
## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpQoXYeb/downloaded_packages
## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpQoXYeb/downloaded_packages
load(here("jk_code", "SO4analysis.rds"))

2 Isolating SO1 and SO4 samples

DimPlot(SO4,split.by = "treatment",group.by = "subclass_MD" )

Idents <- "treatment"
DimPlot(SO4)

Idents(SO4) <- "sample"

DimPlot(SO4,split.by = "sample")

SO <- subset(SO4, idents = c("SO1","SO4"))

DimPlot(SO)

rm(SO4)

3 Finding Markers that define each Control Sample

SO1 <- SCTransform(SO) %>%
    RunPCA() %>%
    FindNeighbors(dims = 1:30) %>%
    FindClusters() %>%
    RunUMAP(dims = 1:30)
## Running SCTransform on assay: RNA
## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## vst.flavor='v2' set. Using model with fixed slope and excluding poisson genes.
## Calculating cell attributes from input UMI matrix: log_umi
## Variance stabilizing transformation of count matrix of size 21397 by 5010
## Model formula is y ~ log_umi
## Get Negative Binomial regression parameters per gene
## Using 2000 genes, 5000 cells
## Found 126 outliers - those will be ignored in fitting/regularization step
## Second step: Get residuals using fitted parameters for 21397 genes
## Computing corrected count matrix for 21397 genes
## Calculating gene attributes
## Wall clock passed: Time difference of 11.76 secs
## Determine variable features
## Centering data matrix
## Place corrected count matrix in counts slot
## Warning: The `slot` argument of `SetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Different cells and/or features from existing assay SCT
## Set default assay to SCT
## PC_ 1 
## Positive:  Pappa2, Malat1, S100g, Zfand5, Pde10a, Ramp3, Robo2, Sdc4, Itga4, Gm37376 
##     Aard, Nadk2, Nos1, Ranbp3l, Wwc2, Sgms2, Neat1, Irx1, Col4a3, Col4a4 
##     Tmem158, Aebp1, Cdkn1c, Itprid2, Cacna1d, Ptgs2, Ccdc80, Bmp3, Camk2d, Uroc1 
## Negative:  Umod, Egf, Tmem52b, Fabp3, Klk1, Sostdc1, Sult1d1, Prdx5, Wfdc15b, Ly6a 
##     Foxq1, Cox6c, Krt7, Ckb, Cldn19, Cox5b, Slc25a5, Wfdc2, Atp5g1, Atp1a1 
##     Ggt1, Cox4i1, Atp1b1, Gm47708, Gadd45g, Cox7b, Ndufa4, Cox8a, Atp5b, Gm44120 
## PC_ 2 
## Positive:  Ubb, Mt1, Fth1, Gm22133, Hspa1a, Ldhb, Actb, Hspa1b, H3f3b, Fos 
##     Junb, Ftl1, Eif1, Jund, Jun, Cd63, Fxyd2, Prdx1, Mt2, Tmem213 
##     Mpc2, Cd9, Mgst1, Rps24, Rpl37, Ier2, Cox8a, Btg2, Gm8355, Egr1 
## Negative:  Mir6236, Gm37376, CT010467.1, Egf, Malat1, Umod, Gm24447, Nme7, mt-Nd5, mt-Rnr1 
##     Slc12a1, mt-Nd6, Kcnq1ot1, Neat1, Tmem52b, Lars2, Atp1b1, mt-Nd4l, Dst, Etnk1 
##     Wnk1, mt-Rnr2, Atrx, Slc5a3, Syne2, Robo2, mt-Co1, mt-Nd2, Zbtb20, Pnisr 
## PC_ 3 
## Positive:  Gm22133, Ldhb, Mgst1, Car15, Mpc2, Cd63, Tmem213, Fxyd2, Cd9, Cldn10 
##     Ftl1, Rpl9, Prdx1, Rplp1, Wfdc2, Rpl7, S100a1, Fth1, Ramp3, Rps27a 
##     Rps24, Ly6a, Irx1, Ppp1r1a, Atp5e, Apoe, Mdh1, Cox8a, Tmbim6, Rpl37 
## Negative:  Fos, Jun, Junb, Egr1, Hspa1b, Btg2, Hspa1a, Atf3, Zfp36, Klf6 
##     Ier2, Fosb, Socs3, Klf2, Jund, Dnajb1, Gadd45g, Tsc22d1, Gm37376, Rhob 
##     Malat1, Gadd45b, Actb, H1f2, Ubc, H2bc4, H3f3b, Egf, Nr4a1, Dusp1 
## PC_ 4 
## Positive:  S100g, Gm8420, Actb, Sdc4, Pth1r, Igfbp7, Gm8885, Pappa2, Malat1, Umod 
##     Ramp3, Tmem52b, Egf, Cebpd, Cd9, Tmsb10, Sgms2, Wnk1, Gm28037, Rpl41 
##     Ppm1h, Tmem158, Ccdc80, Lmna, Tbxas1, Pdcd4, Ppp1r1a, Tmbim6, Col4a3, Gm11808 
## Negative:  Mt1, Mt2, Rpl15-ps2, Rpl13, Ptger3, CT010467.1, Rpl10, Rpl9, Gm22133, mt-Rnr2 
##     Apoe, Rpl24, Aebp1, Rpl17, Jun, Rpsa, Ckb, Hspa8, Eef1a1, Rpl18 
##     Cd63, Gm13340, Gpx6, Gm5905, Uba52-ps, Gm23935, Gpx4, Rplp1, Mfsd4a, Chchd10 
## PC_ 5 
## Positive:  Pappa2, Tmem52b, Umod, Aard, Egf, Sult1d1, Foxq1, Actb, Tmem158, Gm44120 
##     Mcub, Pde10a, Rpl15-ps2, Defb42, Tmsb4x, Dctd, Car15, Rpl10, Kcnj10, 5330417C22Rik 
##     Begain, Rpl13, Rpl24, Abca13, Gsn, Zfand5, Eef1a1, Fabp3, 1700028P14Rik, Gm22133 
## Negative:  mt-Nd5, mt-Co1, mt-Nd4l, Gm28437, mt-Cytb, Malat1, Gm28438, mt-Rnr1, Kng1, Nme7 
##     mt-Nd2, Mir6236, Gm8420, Hspa1b, Gm28037, Gm10925, Defb1, Rps21, Atp1b1, Fxyd2 
##     Klk1, Gm11808, Gm8885, Tnfaip2, Cox6a1, Gm12338, S100g, Sostdc1, mt-Nd3, Car12
## Computing nearest neighbor graph
## Computing SNN
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 5010
## Number of edges: 184583
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7657
## Number of communities: 10
## Elapsed time: 0 seconds
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 11:42:13 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:42:13 Read 5010 rows and found 30 numeric columns
## 11:42:13 Using Annoy for neighbor search, n_neighbors = 30
## 11:42:13 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:42:13 Writing NN index file to temp file /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpQoXYeb/file157107176d16a
## 11:42:13 Searching Annoy index using 1 thread, search_k = 3000
## 11:42:14 Annoy recall = 100%
## 11:42:14 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 11:42:14 Initializing from normalized Laplacian + noise (using RSpectra)
## 11:42:14 Commencing optimization for 500 epochs, with 212904 positive edges
## 11:42:14 Using rng type: pcg
## 11:42:17 Optimization finished
DimPlot(SO1)

DimPlot(SO1,split.by = "sample")

Idents(SO1) <- "sample"
SO1.markers <- FindAllMarkers(SO1, only.pos = TRUE)
## Calculating cluster SO1
## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session
## Calculating cluster SO4
SO1.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)
SO1.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10
DoHeatmap(SO1, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(SO1, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay:
## Ldhb-ps, Gm8292, Rps6

4 Visualizing Difference

FeaturePlot(SO1,"Ptger3")
## Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

FeaturePlot(SO1, "S100g")

VlnPlot(SO1,"Hspa8")

VlnPlot(SO1, "Ptger3")

VlnPlot(SO1, "S100g")

markers.to.plot1 <- c("Ptger3","Hspa8","S100g")
  

DotPlot(SO1,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale = FALSE,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()

SO1 expresses highly in Ptger3, SO4 expresses high in S100g.

S01 Also has a lot more Ribosomal gene influence.

Why?

5 Separating the markers

SO1_markers <- SO1.markers[SO1.markers$cluster == "SO1", ]
SO4_markers <- SO1.markers[SO1.markers$cluster == "SO4", ]

6 View pathway analysis for the genes

6.1 SO1, Control sample 1

df<- SO1_markers %>% arrange(desc(avg_log2FC))

df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% rownames_to_column(var="SYMBOL")

head(markers, n = 50)
ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 4.37% of input gene IDs are fail to map...
markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC)))
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go
## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:1619] "20564" "380785" "94283" "15080" "54377" "622570" "69524" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...1028 enriched terms found
## 'data.frame':    1028 obs. of  12 variables:
##  $ ID            : chr  "GO:0002181" "GO:0022613" "GO:0006457" "GO:0008380" ...
##  $ Description   : chr  "cytoplasmic translation" "ribonucleoprotein complex biogenesis" "protein folding" "RNA splicing" ...
##  $ GeneRatio     : chr  "58/1431" "92/1431" "49/1431" "79/1431" ...
##  $ BgRatio       : chr  "171/28928" "452/28928" "185/28928" "464/28928" ...
##  $ RichFactor    : num  0.339 0.204 0.265 0.17 0.261 ...
##  $ FoldEnrichment: num  6.86 4.11 5.35 3.44 5.27 ...
##  $ zScore        : num  17.5 15.2 13.6 12.1 13.3 ...
##  $ pvalue        : num  6.44e-33 1.27e-31 1.28e-22 4.21e-22 7.06e-22 ...
##  $ p.adjust      : num  3.66e-29 3.63e-28 2.42e-19 5.99e-19 8.03e-19 ...
##  $ qvalue        : num  2.65e-29 2.62e-28 1.75e-19 4.33e-19 5.81e-19 ...
##  $ geneID        : chr  "Gm6133/Rps6/Rpl13/Rpl18/Eif2s3y/Rpl17/Rpsa/Rps19/Rpl9/Rpl10a/Rpl24/Rpl28/Rpl18a/Rpl21/Rpl27a/Rps6-ps4/Rpl12/Rpl"| __truncated__ "Rps6/Eif2s3y/Nol3/Mrto4/Bysl/Rpsa/Nop56/Celf2/Rps19/Rpl24/Utp14b/Fbl/Zfp593/Gemin7/Mphosph10/Rps6-ps4/Srsf6/Fts"| __truncated__ "Hspa8/Fkbp11/Ric3/Umod/Hspb1/Hsph1/Ranbp2/Dnajc2/Dnajc10/Ahsa2/Tor2a/Nktr/Ppil1/Mkks/Hspa4/Chordc1/Cct3/Aip/Dna"| __truncated__ "Hspa8/Nol3/Celf2/Tsen34/U2af1/Umod/Srsf3/Polr2a/Gemin7/Srsf6/Srsf7/Zc3h13/Ybx1/Mfap1a/Prpf8/Zfp326/Rnpc3/Phf5a/"| __truncated__ ...
##  $ Count         : int  58 92 49 79 48 63 48 74 64 23 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 lines
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

I Think that SO1, has too much ribosomal influence so I will try to remove it here and see what the changes are.

6.2 SO4 Control Sample 2

df<- SO4_markers %>% arrange(desc(avg_log2FC))

df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% rownames_to_column(var="SYMBOL")

head(markers, n = 50)
ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 8.82% of input gene IDs are fail to map...
markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC)))
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go
## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:424] "667014" "100041678" "667937" "671215" "13653" "384179" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...377 enriched terms found
## 'data.frame':    377 obs. of  12 variables:
##  $ ID            : chr  "GO:0140236" "GO:0140241" "GO:0140242" "GO:0002181" ...
##  $ Description   : chr  "translation at presynapse" "translation at synapse" "translation at postsynapse" "cytoplasmic translation" ...
##  $ GeneRatio     : chr  "14/306" "14/306" "14/306" "21/306" ...
##  $ BgRatio       : chr  "47/28928" "48/28928" "48/28928" "171/28928" ...
##  $ RichFactor    : num  0.2979 0.2917 0.2917 0.1228 0.0806 ...
##  $ FoldEnrichment: num  28.16 27.57 27.57 11.61 7.62 ...
##  $ zScore        : num  19.27 19.05 19.05 14.39 9.37 ...
##  $ pvalue        : num  4.07e-17 5.69e-17 5.69e-17 1.67e-16 1.54e-09 ...
##  $ p.adjust      : num  7.33e-14 7.33e-14 7.33e-14 1.61e-13 1.19e-06 ...
##  $ qvalue        : num  5.46e-14 5.46e-14 5.46e-14 1.20e-13 8.85e-07 ...
##  $ geneID        : chr  "Rpl38/Rps28/Rpl37a/Rpl37/Rps26/Rps27/Rpl23a/Rpl36a/Rpl34/Rps27a/Rpl36/Rps12/Rpl22/Rpl35a" "Rpl38/Rps28/Rpl37a/Rpl37/Rps26/Rps27/Rpl23a/Rpl36a/Rpl34/Rps27a/Rpl36/Rps12/Rpl22/Rpl35a" "Rpl38/Rps28/Rpl37a/Rpl37/Rps26/Rps27/Rpl23a/Rpl36a/Rpl34/Rps27a/Rpl36/Rps12/Rpl22/Rpl35a" "Rpl38/Rps28/Rpl37a/Rpl37/Rpl41/Rpl39/Rps29/Rps21/Rpl31/Pabpc1/Rps26/Rpl22l1/Rpl6l/Rpl23a/Rpl36a/Rpl34/Rps27a/Rp"| __truncated__ ...
##  $ Count         : int  14 14 14 21 15 16 21 19 18 19 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 lines
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.